#######################################
######Masculinity Threat Analysis######
#######################################
###############10/10/25##################
#######################################

####Analysis Script for Gothreau & Haas, forthcoming in Journal of Experimental Political Science###
####A Replication and Extension of Willer et al. (2013) Overdoing Gender: A Test of the Masculine Overcompensation Thesis### 

library(haven)
library(ggplot2)
library(dplyr)
library(broom)
library(officer)
library(flextable)
library(forcats)
library(psych)
library(survey)
library(tidyr)
library(stargazer)
library(forcats)

#Set Working Directory#
MTData <- read.csv("CleanMTData.csv")

#####EXPLORATORY FACTOR ANALYSIS#####

fa_vars <- c(
  "IRAQ_WAR_rescaled", "GAY_RIGHT_rescaled", "SUV_DESIRABILITY_rescaled", 
  "IMMIGRATION_rescaled", "ELECTRIC_CAR_rescaled", "TRANS_RIGHT_rescaled", 
  "SUV_PAY_rescaled", "HIRINGPOLICY_WOMEN_rescaled", "SYSTEM_JUSTIFICATION_rescaled",
  "MARIJUANA_rescaled", "CONSERVATISM_rescaled", 
  "DOMINANCE_scale", "TRADITIONALISM_scale"
)
fa_data <- MTData %>%
  select(all_of(fa_vars)) %>%
  na.omit()  

#KMO measure of sampling adequacy 
KMO(fa_data)  # >0.6 
cortest.bartlett(fa_data)

#Scree plot to decide how many factors to extract 
fa.parallel(fa_data, fm = "ml", fa = "fa")  # 'ml' = maximum likelihood method

efa_result <- fa(fa_data, nfactors = 2, rotate = "oblimin", fm = "ml")
print(efa_result, sort = TRUE)

fa.diagram(efa_result)

#Clean labels
item_labels <- c(
  IRAQ_WAR_rescaled             = "Iraq War",
  GAY_RIGHT_rescaled            = "Gay Rights",
  SUV_DESIRABILITY_rescaled     = "SUV Desirability",
  IMMIGRATION_rescaled          = "Immigration",
  ELECTRIC_CAR_rescaled         = "Electric Car",
  TRANS_RIGHT_rescaled          = "Trans Rights",
  SUV_PAY_rescaled              = "SUV Payment",
  HIRINGPOLICY_WOMEN_rescaled   = "Gender Affirmative Action",
  SYSTEM_JUSTIFICATION_rescaled = "System Justification",
  MARIJUANA_rescaled            = "Marijuana Legalization",
  CONSERVATISM_rescaled         = "Conservatism",
  DOMINANCE_scale               = "Social Dominance",
  TRADITIONALISM_scale          = "Traditionalism"
)

#Loadings into a data frame
loadings_df <- as.data.frame(unclass(efa_result$loadings))
loadings_df$Item <- item_labels[rownames(loadings_df)]
loadings_df <- loadings_df %>%
  select(Item, everything()) %>%
  rename(`Factor 1` = ML1, `Factor 2` = ML2)

#Round loadings
ft_df <- loadings_df %>%
  mutate(`Factor 1` = round(`Factor 1`, 2),
         `Factor 2` = round(`Factor 2`, 2))

#Create flextable
ft <- flextable(ft_df)
ft <- set_caption(ft, "Factor Loadings from Exploratory Factor Analysis")
ft <- autofit(ft)

ft <- bold(ft, i = ~ `Factor 1` >= 0.4, j = "Factor 1", bold = TRUE)
ft <- bold(ft, i = ~ `Factor 2` >= 0.4, j = "Factor 2", bold = TRUE)

doc <- read_docx()
doc <- body_add_par(doc, "Table X. Factor Loadings from Exploratory Factor Analysis", style = "heading 1")
doc <- body_add_flextable(doc, ft)
print(doc, target = "factor_loadings_table.docx")

#######Generate factor scores 

MTData_fa <- MTData %>%
  filter(complete.cases(select(., all_of(fa_vars))))

factor_scores <- factor.scores(MTData_fa[, fa_vars], efa_result, method = "regression")

MTData_fa$Factor1 <- factor_scores$scores[, 1]
MTData_fa$Factor2 <- factor_scores$scores[, 2]

#split by gender
Data_male <- MTData_fa %>% filter(gender_binary == "Men")
Data_female <- MTData_fa %>% filter(gender_binary == "Women")

####T-Tests

#Define DVs and conditions
factor_scores <- c("Factor1", "Factor2")
control_group <- "no threat"
other_conditions <- c("threat", "constrained", "PCKT")

#Function to run t-tests
run_gender_ttests <- function(data, gender_label) {
  results <- list()
  for (score in factor_scores) {
    for (cond in other_conditions) {
      test <- t.test(data[[score]][data$P_EXP == cond],
                     data[[score]][data$P_EXP == control_group])
      results[[paste(gender_label, score, cond, sep = "_vs_")]] <- tidy(test)
    }
  }
  return(results)
}

ttests_male <- run_gender_ttests(Data_male, "Men")
ttests_female <- run_gender_ttests(Data_female, "Women")

ttests_combined <- do.call(rbind, c(ttests_female, ttests_male))
ttests_combined$comparison <- rownames(ttests_combined)
rownames(ttests_combined) <- NULL

####Create Table (Table C4 in the Appendix)

ttests_clean <- separate(ttests_combined, comparison, into = c("gender", "DV", "experimental_group"), sep = "_vs_")
ttests_clean <- ttests_clean %>%
  select(gender, DV, experimental_group, estimate, statistic, p.value, conf.low, conf.high) %>%
  rename(
    t = statistic,
    p = p.value,
    CI_low = conf.low,
    CI_high = conf.high
  ) %>%
  mutate(across(where(is.numeric), round, 3))

tt_ft <- flextable(ttests_clean)
tt_ft <- autofit(tt_ft)
doc <- read_docx()
doc <- body_add_par(doc, "T-Tests for Factor Scores by Gender and Experimental Condition", style = "heading 1")
doc <- body_add_flextable(doc, tt_ft)
print(doc, target = "factor_score_ttests.docx")

#####Regression Models 
#controls for respondent gender, age, education, and ideology 

#Create treatment indicator variables 
MTData <- MTData %>%
  mutate(
    P_EXP = factor(P_EXP),                  
    P_EXP = fct_relevel(P_EXP, "no threat")        
  )

Data_female <- subset(MTData, gender_binary == "Women")
Data_male   <- subset(MTData, gender_binary == "Men")

MTData$P_EXP <- relevel(MTData$P_EXP, ref = "no threat")
Data_female$P_EXP <- relevel(Data_female$P_EXP, ref = "no threat")
Data_male$P_EXP<- relevel(Data_male$P_EXP, ref = "no threat")

#Regression Models with interactions 

M1 <- lm(IRAQ_WAR_rescaled ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M2 <- lm(GAY_RIGHT_rescaled ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M3 <- lm(SUV_DESIRABILITY_rescaled ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M4 <- lm(SUV_PAY_rescaled ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M5 <- lm(DOMINANCE_scale ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M6 <- lm(SYSTEM_JUSTIFICATION_rescaled  ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M7 <- lm(CONSERVATISM_rescaled ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M8 <- lm(IMMIGRATION_rescaled  ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M9 <- lm(MARIJUANA_rescaled  ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M10 <- lm(ELECTRIC_CAR_rescaled  ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M11 <- lm(TRANS_RIGHT_rescaled ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)
M12 <- lm(HIRINGPOLICY_WOMEN_rescaled ~ P_EXP + P_EXP*gender_binary + gender_binary + AGE + IDEO + EDUC5, data=MTData, weights = WEIGHT)

#Export to Word document 
stargazer(M1, M2, M3, M4, M5, M6, M7, M8, M9, M10, M11, M12,
          type = "html",                      
          out = "regression_results.html",    
          title = "Regression Results by Outcome Variable",
          dep.var.labels = c("Iraq War", "Gay Rights", "SUV Desirability", "SUV Payment",
                             "Dominance", "System Justification", "Conservatism",
                             "Immigration", "Marijuana", "Electric Car",
                             "Trans Rights", "Hiring Policy (Women)"),
          covariate.labels = c("Threat", "Constrained", "PCKT",
                               "Female",
                               "Threat x Female", "Constrained x Female", "PCKT x Female",
                               "Age", "Ideology", "Education"),
          omit.stat = c("f", "ser"),  
          align = TRUE,
          no.space = TRUE,
          single.row = TRUE)


#FEMALE models
M1_f <- lm(IRAQ_WAR_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M2_f <- lm(GAY_RIGHT_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M3_f <- lm(SUV_DESIRABILITY_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M4_f <- lm(SUV_PAY_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M5_f <- lm(DOMINANCE_scale ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M6_f <- lm(SYSTEM_JUSTIFICATION_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M7_f <- lm(CONSERVATISM_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M8_f <- lm(IMMIGRATION_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M9_f <- lm(MARIJUANA_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M10_f <- lm(ELECTRIC_CAR_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M11_f <- lm(TRANS_RIGHT_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)
M12_f <- lm(HIRINGPOLICY_WOMEN_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_female, weights = WEIGHT)

#MALE models
M1_m <- lm(IRAQ_WAR_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M2_m <- lm(GAY_RIGHT_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M3_m <- lm(SUV_DESIRABILITY_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M4_m <- lm(SUV_PAY_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M5_m <- lm(DOMINANCE_scale ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M6_m <- lm(SYSTEM_JUSTIFICATION_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M7_m <- lm(CONSERVATISM_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M8_m <- lm(IMMIGRATION_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M9_m <- lm(MARIJUANA_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M10_m <- lm(ELECTRIC_CAR_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M11_m <- lm(TRANS_RIGHT_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)
M12_m <- lm(HIRINGPOLICY_WOMEN_rescaled ~ P_EXP + AGE + IDEO + EDUC5, data = Data_male, weights = WEIGHT)

#Export to Word document 
# === FEMALE TABLE ===
stargazer(M1_f, M2_f, M3_f, M4_f, M5_f, M6_f,
          type = "html",
          out = "regression_results_female.html",
          title = "Regression Results: Female Respondents")

stargazer(M7_f, M8_f, M9_f, M10_f, M11_f, M12_f,
          type = "html",
          out = "regression_results_female_2.html",
          title = "Regression Results: Female Respondents")

stargazer(M1_m, M2_m, M3_m, M4_m, M5_m, M6_m,
          type = "html",
          out = "regression_results_male.html",
          title = "Regression Results: Male Respondents")

stargazer(M7_m, M8_m, M9_m, M10_m, M11_m, M12_m,
          type = "html",
          out = "regression_results_male_2.html",
          title = "Regression Results: Male Respondents")



